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1. Introduction 

The observation of parity-time (VT) symmetry in experiments with structured optical 
waveguides [TJE] marked a remarkable triumph for the mathematical theory of non- 
Hcrmitian Hamiltonians having VT symmetry, initiated by Bender and Bocttcher [3J 
in 1999. The experiments exploit the optical-quantum mechanical analogy that the 
wave equation for the transverse electric field mode is formally equivalent to the 
one-dimensional Schrodinger equation when the potential is identified with V(x) = 
— k 2 n 2 (x) and the energy eigenvalue with E = —/3 2 , where k = u>/c is the vacuum wave 
number and /3 the propagation constant of the mode. A VT symmetric experiment is 
realized by introducing loss and gain terms in the complex index of refraction n(x). 

A quantum analogue to the optical experiments would be a Bose-Einstein 
condensate in a double-well potential where loss and gain is realized by removing 
atoms from one well and pumping in atoms in the other, as has been suggested by 
Klaiman et al. [4] . A complication arises from the interaction between the condensate 
atoms. In the simplest case this is s-wave scattering, which yields an additional (real- 
valued) potential term proportional to \ip\ 2 , and turns the non-Hermitian Hamiltonian 
into a nonlinear one. For the Hamiltonian to be VT symmetric this term has to be a 
symmetric function of the coordinates, which cannot be guaranteed from the outset. 
It will be the purpose of this paper to analyze whether or not the nonlinearity destroys 
the VT symmetry of the quantum system. 

We will do this, in the spirit of a model calculation, by considering the situation 
where the double well is idealized by two attractive delta-function potential wells, with 
loss added in one trap and gain in the other. We will demonstrate that the stationary 
solutions of the Gross-Pitaevskii equation indeed preserve the VT symmetry of the 
nonlinear Hamiltonian, and merge in a branch point at some critical value of the loss 
and gain, beyond which the VT symmetry is broken. 

Another model which simulates the physical situation of a Bose-Einstein 
condensate in a double well with loss and gain has been investigated by Graefe et al. 
[5HZ] in the framework of a two-mode Bosc-Hubbard-type VT symmetric Hamiltonian. 
In the model the dynamics of the quantum system is mapped on that of a (formal) spin 
vector whose motion on the surface of the Bloch sphere can be analyzed in terms of 
classical nonlinear dynamics. As an optical analogue, in the two-mode approximation 
Ramezani et al. [8] have recently looked at a mathematical model of a VT symmetric 
coupled dual waveguide arrangement with Kerr nonlinearity. It will be interesting to 
see which features of these models are recovered when actually solving the nonlinear 
VT symmetric Gross-Pitaevskii equation. 

Delta-function potentials are popular as model systems in the literature since 
they allow for analytic or partially analytic solutions but are flexible enough to 
provide insight into characteristic phenomena of more complex physical situations. 
For example, resonances and the decay behaviour of a Bose-Einstein condensate in a 
double delta-shell potential were analyzed by Rapedius and Korsch [9] with the aim of 
modelling features typical of double- well potentials [TU] ■ In the same vein, Witthaut et 
al. have studied the nonlinear Schrodinger equation for a delta-functions comb to 
gain insight into the properties of nonlinear stationary states of periodic potentials. For 
a system where a real delta-function potential is augmented by a VT symmetric pair of 
delta-functions with imaginary coefficients, bound states and scattering wave functions 
have been calculated by Jones [12| . The paper was devoted to the quasi- Hermitian 
analysis of the problem, and no nonlinearity was present. Jakubsky and Znojil |13j 
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have considered the explicitly solvable model of a particle exposed to two imaginary 
VT delta-function potentials in an infinitely high square well, and determined the 
energy spectrum, but also without nonlinearity. A large body of literature of course 
exists on the properties of solitons and Bose-Einstein condensates in periodic optical 
and nonlinear lattices with VT symmetry and their nonlinear optical analogues (see, 
e.g., [T4TI22"] ) . But to the best of our knowledge the basic problem of the nonlinear 
Schrodinger equation with two VT symmetric delta-function double wells has not been 
treated so far. 

2. Delta-function traps model 

We consider a Bose-Einstein condensate trapped in two delta-function shaped 
potential wells located at x = ±a/2, where from one well condensate atoms are 
removed and the same amount of atoms is added in the other. At sufficiently low 
temperatures stationary solutions for the condensate wave function are then described 
in a mean-field approach by the nonlinear Schrodinger equation, or Gross-Pitaevskii 
equation (see, e.g., [23]). 

- V"(x) - [(1 + i-y)5{x + a/2) + (1 - i 7 )<5(x - a/2)] #(x) 

-g\y(x)\ 2 y(x) = ~K 2 y(x), (1) 

where g is the nonlinear interaction strength, and the real-valued parameter 7 
determines the size of the gain and loss terms. Units have been chosen in such a 
way that the strength of the real part of the delta-function potentials is normalized to 
unity. (The relation to physical quantities is given in the Appendix.) For stationary 
solutions the eigenvalues n will be real, but since we are interested in the complete 
eigenvalue spectrum, including complex eigenvalues, we will quite generally search 
for solutions with k e C, Re(n) > 0. While the delta-function potentials are VT 
symmetric, it is not clear a priori that the equation itself is VT symmetric since this 
requires the nonlinear potential term to be a symmetric function of x. 

To be in a position to study the effects that emerge when the nonlinearity is turned 
on we first consider the case where the nonlinearity is absent, and then proceed to 
values g > 0. 

2.1. Delta- functions trap model for vanishing nonlinearity 

For 3 = the bound-state wave function has the form: 

!Ae KX : x < -b 

Ce KX + De- KX : -b < x < b 
Be~ KX : b < x 

Applying the continuity conditions leads to the system of linear equations (with the 
abbreviation kq = 1 + ij) , 

n e- Kb C+(k q - 2n)c Kb D = , (2) 
(K*-2K)e Kb C + K*e- Kb D = 0. (3) 

The eigenvalues k are obtained numerically by finding the roots of the complex secular 
equation 

(l + 7 2 )c- 2Ka -(l+ 7 2 +4«; 2 -4K)=0. (4) 
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Figure 1. Eigenvalues re of equation |(T} for vanishing nonlinearity as functions of 
the size of the loss/gain parameter 7 for different distances of the delta-function 
traps: (a) a = 2.2, (b) a = 2.0. (c) a = 1.8. 



It depends on two parameters, the distance a of the traps and the strength 7 of the 
loss/gain terms. 

In figure Q] the eigenvalues k are plotted as functions of the strength of the 
loss/gain parameter 7 for different distances of the traps. It can be seen that the 
solutions exhibit the behaviour typical of VT symmetric non-Hermitian Hamiltonians: 
up to a critical value of the loss and gain parameter 7 there exist two real eigenvalues 
K. Since we have included a minus sign on the right-hand side of the eigenvalue 
equation ([T]), in the range of real eigenvalues the upper branch belongs to the 
ground state and the lower branch to an excited state. At the critical value - 
the branch point (or exceptional point) - the eigenvalues and the eigenfunctions 
merge. Beyond the branch point the eigenvalues become complex, and complex 
conjugate. The time-dependence of the eigenstates of ((T|) with g = is given by 
exp(ift 2 t) = exp(i(fc 2 — kf)t) exp(— 2fc r fcit) (with fc r and k\ the real and imaginary part 
of k), hence the mode with fc; > decays, and the mode with fc, < grows. 

Of course at 7 = the two states turn into the well-known cosh(ra) (ground 
state) and sinh(Kx) (excited state) textbook solutions of ([T]) between the traps. It is 
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Figure 2. Real parts, (a), and imaginary parts, (b), of the ratio of the coefficients 
of the wave functions in the region between the delta-function traps for vanishing 
nonlinearity and a distance of a = 2.2 as functions of the size of the loss/gain 
parameter 7. For values of 7 below the bifurcation point the solid lines correspond 
to the ground state, the dashed lines to the excited state. 



also well-known that for 7 = the excited state disappears at a distance of a — 2. 
From figure [T] (b) we see, however, that for a — 2 an excited state is born as soon as 
loss and gain is switched on. The Figure also shows that for distances below a = 2 
(cf. figure [1] (c)) the excited state is born when loss and gain exceed a finite threshold 
value. 

The behaviour of the eigenvalues shown in figure [T] is in complete analogy with 
that found by Klaiman et al. in their study of two micro waveguides with loss and 
gain (cf. figure 2 in Ref. [4]). Their calculation was performed for a fixed distance 
of the waveguides. It would be interesting to investigate whether or not the excited 
mode also vanishes when the distance of the waveguides is decreased. 

Below the branch point the wave functions possess VT symmetry, *£*(x) = 
^(—x). In our model this can easily be seen since for real k the complex conjugate 
of © is identical to ^ so that the pairs (C,D) and (D*,C*) fulfill the same 
equation, and, hence D = C* . From the continuity condition then follows that 
A = C + C*e 2Kb and B = Ce 2Kb + C* = A*. Note that quite generally the VT 
symmetry of a wave function implies that its modulus is a symmetric function, since 
= *(-x)**(-af). 

Figure [2] shows, for a = 2.2, the ratio of the coefficients C and D of the wave 
functions in dependence of the loss/gain parameter 7. What is shown is the evolution 
of the modulus of the ratio and its phase for the ground state and the excited state. 
Since below the branch point we have D — C* — |C|exp(— if), the modulus of the 
ratio is 1, and the phase is twice the phase ip of the coefficient C. As it must be, at 
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Figure 3. Time evolution of a superposition of the two eigenstates for different 
values of the loss/time parameter: (a) 7 = 0: (b) 7 = 0.35, (c) 7 = 0.39, close to 
the exceptional point. 



7 = the phase starts with a value of for the (symmetric) ground state and it for the 
(antisymmetric) excited state. While the phase for the ground state is quickly rotated 
away from its initial value, the phase of the ratio belonging to the excited state only 
slowly deviates from ir up to the branch point. At the branch point the modulus and 
the phase of D/C coincide. Beyond the branch point they quickly evolve away from 
the values they have in the VT symmetric regime. 

It is instructive to look at the time dependence of a superposition of the two 
eigenstates in the VT symmetric regime. For different values of the loss/gain 
parameter hgure [3] shows the time evolution of 

\^(x,t)\ 2 = \i>\{x) exp(i«?i) + tp 2 (x) exp(inlt)\ 2 . (5) 

Without the loss/gain term the probability density oscillates between the two delta- 
function wells. For finite gain and loss the beat frequency decreases and the oscillation 
is deformed, and for a loss/gain parameter close to the branch point the probability 
density on both sides is almost identical and pulsates. This behaviour completely 
agrees with the time evolution of the power distribution for a propagating electric 
sum field consisting of two guided modes calculated by Klaiman et al. (cf. figure 3 in 
Ref. 

Thus we find that the simple quantum mechanical model studied in this section 
already captures, for both eigenvalues and wave functions, all the effects of a VT 
symmetric wave guide configuration in optics. 
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2.2. Delta- functions trap model with nonlinearity 

2.2.1. Numerical method For j^Owe have solved the Gross-Pitaevskii equation (JXJ) 
numerically using a procedure in which the energy eigenvalues are found by a five- 
dimensional numerical root search. The free parameters which have to be adjusted in 
such a way that a physically meaningful wave function is obtained are the eigenvalue k 
as well as initial conditions for the wave function and its derivative. Since the overall 
phase is arbitrary we can choose it such that ^(O) is a real number. Therefore five 
real parameters remain, viz. the real part of ^(0), and the real and imaginary parts 
of both v]>'(0) and n. Physically relevant wave functions must be square integrable 
and normalized. The normalization is important since the Gross-Pitaevskii equation 
is nonlinear and the norm influences the form of the Hamiltonian. This gives in total 
five conditions which have to be fulfilled: The real and imaginary parts of ^ must 
vanish for x — > ±oo, and the norm of the wave function must fulfill — 1 = 0. 

Outside the delta-function traps the Gross-Pitaevskii equation (TTJ) coincides with 
the free nonlinear Schrodinger equation, which has well known real solutions in terms of 
Jacobi elliptic functions (cf., e.g., [11,24,25 ). The function which solves the equation 
in the ranges |x| > a/2 for the attractive nonlinearity considered here and decays to 
zero for |x| — > oo is cn(/tx, 1) = sech(Kx) = 1/cos1i(kx). We find that once the correct 
eigenvalues and eigenfunctions are obtained our numerical wave functions exactly show 
this behaviour. 

2.2.2. Eigenvalues Figure [4] shows the results for the real and imaginary parts of 
the eigenvalues k calculated for a trap distance of a = 2.2 and different values of the 
nonlinearity parameter as functions of 7. The results for the case g = are shown for 
comparison. It can be recognized that even with the nonlinearity present there always 
exist two branches of real eigenvalues, up to critical values 7 cr , where the two branches 
merge. It can also be seen that branches of two complex conjugate eigenvalues also 
appear, but it comes as a surprise that these are born, not at 7 cr , but at smaller 
values of 7 where they bifurcate from the real eigenvalue branch of the ground state. 
As the nonlinearity is increased, the bifurcation points are shifted to smaller values 
of 7, while the branches of the real eigenvalues practically retain their form and are 
shifted upwards. We therefore find that for nonvanishing nonlinearity there always 
exists a range of 7 values where two real and two complex eigenvalues coexist. 

This behaviour obviously is a new effect, caused by the nonlinearity in the Gross- 
Pitaevskii equation. In fact the different branches can be continuously transformed 
into each other: when, for fixed 7, we increase the nonlinearity from g = in small 
steps to g — 1 and always take the solution of the previous step as input for the root 
search in the next step, the branch of complex conjugate eigenvalues which exists for 
g = is continuously transferred to the branch of complex eigenvalues for g — 1 shown 
in figure 01 Likewise, starting from the branch of complex conjugate eigenvalues for 
g = 1 at any value of 7 in the range where four solutions exist and decreasing g in 
small steps to g = we end up on the branch belonging to the real eigenvalue of the 
ground state. 

At this point it is useful to compare with the results of the investigations of a 
VT symmetric two-mode Bose-Hubbard Hamiltonian with loss and gain by Graefe et 
al. [7] . An eigenenergy spectrum with a structure similar to the one shown in figure |4] 
also appeared in their calculations (see figure 13 in 7 ). In the model, stationary states 
correspond to fixed points of the motion of a vector on the Bloch sphere, whose types 




Figure 4. Real parts (a) and imaginary parts (unless zero) (b) of the eigenvalues 
k of the full nonlinear equation {T} as functions of the loss/gain parameter 7 for 
a = 2.2 for the four values of the nonlinearity parameter (in (a) from bottom to 
top) g = 0.0, 0.1, 0.5, and 1.0. Solid lines denote purely real eigenvalues, dashed 
lines complex conjugate eigenvalues. For g ^ the complex conjugate eigenvalues 
bifurcate from the ground state branch before the branch point where the two real 
solutions coalesce. 



can be classified according to the eigenvalues of the Jacobian matrix. In the region 
where only two real eigenvalues exist the fixed points can be identified as centres, 
while in the region with four eigenvalues the fixed points correspond to a centre and 
a saddle point, and a sink and a source. The centre and saddle point collide at the 
branch point and vanish. 

This behaviour is in beautiful agreement with the results shown in figure [U 
A stability analysis of the wave functions shows that up to the bifurcation point 
both the ground state and the excited state are stable, they correspond to the two 
centres. Beyond the bifurcation point the excited state remains stable (a centre) 
while the ground state becomes unstable (a saddle point). Out of the wave functions 
belonging to the two complex conjugate eigenvalues one decays and the other grows, 
corresponding to the sink and the source, respectively. It is gratifying that the two so 
different approaches qualitatively yield the same behaviour. It can be concluded that 
quite generally the familiar branching scheme for linear VT symmetric Hamiltonians 
will be changed into a scheme of the type shown in figure 0] if a nonlinearity is added 
to the Hamiltonian. 



2.2.3. Eigenfunctions In figure [5] we present examples of wave functions of the 
eigenstates of the nonlinear Schrodinger equation ([T|) determined numerically for 
g = 0.5 and a = 2.2. The results for the ground state and the excited state are shown 
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Figure 5. Real and imaginary parts and moduli of wave functions of the 
eigenstates of the nonlinear Hamiltonian in {T} for g = 0.5, a = 2.2. (a) Ground 
state and (b) excited state for 7 = 0.35, (c) solution with imaginary part of k > 
and (d) imaginary part of k < for 7 = 0.5. In (a) and (b), the wave functions 
are VT symmetric, and the moduli are symmetric functions, producing the VT 
symmetry of the total nonlinear Hamiltonian. The wave functions in (c) and (d) 
are no longer VT symmetric, their moduli are not symmetric functions, and the 
VT symmetry also of the nonlinear Hamiltonian is broken. 



in figure[5](a) and (b) for 7 = 0.35, below the the critical value 7 cr « 0.4. The fact that 
in the numerical calculation the global phase of the complex wave function was fixed in 
such a way that V&(0) is real makes it particularly easy to recognize the VT symmetry 
of the two solutions, since their real parts become even functions and their imaginary 
parts odd functions of x. From the VT symmetry of the wave function follows that the 
modulus, also shown in figure [SJ is an even function, and with it the nonlinear term in 
{1}. We therefore have the important result that the nonlinear Schrodinger equation 
with loss and gain selects as eigenfunctions exactly those states in Hilbert space which 
render the nonlinear Hamiltonian VT symmetric! In the ground state, which emerges 
from the symmetric real wave function for 7 = g = 0, the symmetric contribution 
from the real part still dominates, while for the excited state, which originates from 
the antisymmetric solution for 7 = 9 = 0, the antisymmetric contribution from the 
imaginary part prevails. 

The VT symmetry of the wave functions is broken for the eigenstates with 
complex eigenvalues. Figures [5] (c) and (d) show as examples the wave functions 
obtained for 7 = 0.5 for the corresponding pair of complex conjugate eigenvalues k. It 
can be seen that the real and imaginary parts are no longer even or odd functions, and 
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Figure 6. Wave functions for a = 2.2 and g = 0.5 and values of 7 close to 
the exceptional point: ground state (a) and excited state (b) for 7 = 0.395, 
ground state (c) and excited state (d) for 7 = 0.398. While for 7 = 0.395 still 
small differences can be recognized, in particular in the imaginary part, the wave 
functions practically conincide at 7 = 0.398. 



therefore VT symmetry is lost. As a consequence, the moduli of the wave functions 
also are no longer even functions of x. Thus we find that beyond the branch point not 
only the VT symmetry of the wave functions is broken but also that of the nonlinear 
Hamiltonian! 

It must be noted, however, that eigenstates with complex eigenvalues, in spite of 
solving the nonlinear eigenvalue problem ([1]), do not obey the corresponding time- 
dependent nonlinear Schrodinger equation. For complex eigenvalues the modulus 
squared of the wave functions grows or decays proportional to exp(— 2Im(«; 2 )i), and 
so does the nonlinear term in (JT]). Therefore strictly speaking multiplying these 
eigenstates with the usual time evolution factor exp(?K 2 i) only captures the onset 
of the temporal evolution of the two modes, for times Im(/t 2 )£ <C 1. If one wished 
to obtain the exact time evolution of these solutions in the region beyond the branch 
point one would have to solve the time-dependent nonlinear Schrodinger equation (see 
subsection l2.2.5|) . In figure Etc) the mode with positive imaginary part of k is the one 
which decays, as expected it is more strongly localized in the trap with loss, while the 
mode with negative imaginary part in figure [SJd) is the one which grows and is more 
strongly localized in the trap with gain. 

To demonstrate that at 7 cr the two eigenstates with real eigenvalues indeed 
coalesce, we have plotted their wave functions in figure |6] for a = 2.2, g = 0.5, 
close to the critical value of j CI ~ 0.4. While slightly below the critical value minor 
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Figure 7. Eigenvalues of the nonlinear Hamiltonian with a = 2.2 for g = 0.5 
including the analytical continuation as functions of the loss/gain parameter: (a) 
real part, (b) imaginary part of n (unless zero). Now a pair of complex conjugate 
eigenvalues also emerges at the exceptional point. 
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Figure 8. Wave functions of the eigenstates of the nonlinear Hamiltonian 
belonging to the complex eigenvalues of the analytical continuation for g = 0.5 
and a = 2.2 and a loss/gain parameter of 7 = 0.5. The eigenfunctions are 
not VT symmetric. The black dotted lines have been drawn to illustrate that 
lm(*(0)) + 0. 



deviations of the wave functions are still visible, in particular in the imaginary parts, no 
differences can be recognized anymore at 7 = 0.398. Thus also for strong nonlinearity 
at 7 cr we find the properties characteristic of an exceptional point, i.e. the coalescence 
of both eigenvalues and eigenfunctions. 
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2.2.4- Analytical continuation of the nonlinear Hamiltonian The fact that at the 
branch point two real solutions coalesce without giving rise to two solutions with 
complex eigenvalues contradicts the behaviour typical of exceptional points. Obviously 
these solutions cannot be found by solving the nonlinear Gross-Pitaevskii equation 
in its form ([1]), but require an analytical continuation of the nonlinear Hamiltonian 
beyond the critical point 7 cr - The reason is that the nonlinear term pl^l 2 is a 
nonanalytic function, and some care has to be taken when analytically continuing 
the Hamilitonian beyond the exceptional point. 

In the VT symmetric regime up to 7 cr we have = \&(— x). Therefore on the 

way to the bifurcation point we can write the nonlinearity for the VT symmetric states 
in the form t/|\I/(a:) | 2 = g^(x)^(— x). This function can be continued analytically. 
In the numerical calculation the additional condition J ^(x)^(—x)dx = 1 must be 
introduced to fix the phase of the nonlinearity in the VT broken regime. In the 
(then) six-dimensional root search also Iu^^O)) must be varied. As a result we find 
two more complex conjugate solutions that emerge from the coaelescing states, see 
figure [Jj The real and imaginary parts of the eigenfunctions and their moduli are 
shown in figure [H Obviously the states are not VT symmetric, and no longer possess 
vanishing imaginary part at the origin. 

2.2.5. Temporal evolution We finally take a look at the time evolution of a 
superposition of the two stationary eigenstates in the VT symmetric regime. While the 
time dependence of each individual state is given by exp(iK 2 t), the time evolution of 
a superposition can no longer be obtained by simply multiplying the individual states 
by their time evolution factors. The reason is again that the nonlinear term in the 
Gross-Pitaevskii equation becomes time-dependent, and therefore to determine the 
temporal evolution of an initial superposition the time-dependent Gross-Pitaevskii 
equation has to be solved. We have done this numerically by applying the split- 
operator technique, where in the short-time approximation the time evolution operator 
can be split into products of exponentials of the kinetic and potential terms of the 
(nonlinear, time-dependent) Hamiltonian, and the time propagation is calculated by 
successively switching between position and momentum space. 

In figure [S] we show the results for the time evolution of |\&(x, f)| 2 of the initial 
superposition of the ground and the excited state ^(x,0) = (ipi(x) + ^(x)) /y/2, for 
the nonlinearity g = 0.5 and different values of the loss/gain parameter. It can be seen 
that, even with the nonlinearity present, for small 7 we have the behaviour typical 
of VT symmetric systems, i.e. the decrease of the beat length as the loss and gain 
parameter is growing and the appearance of time intervals in which the probability 
density in both wells is simultaneously large, or almost completely extinguished. 

In figure[9](d), for the value of 7 = 0.3 the VT symmetry obviously is broken: even 
though we start with a high probability density in the well with loss, as time proceeds 
there is no beating, and the well with gain wins. Note that this value of 7 is smaller 
than the value of 7 ~ 0.38 where the two complex eigenvalues bifurcate from the 
ground state branch and smaller than the critical value 7 cr 0.4 where the two states 
coalesce. The reason is that the norm of the state constantly varies with time, and 
with it the size of the nonlinear term in the Gross-Pitaevskii equation. This effectively 
corresponds to a variation of g in the stationary Gross-Pitaevskii equation, and acts as 
if one would constantly be moving through different eigenvalue diagrams. Therefore 
for a given value of g, the transition to the VT broken phase of the condensate can 
set in at values of 7 smaller than those of the bifurcation and the exceptional point 
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Figure 9. Temporal evolution of a superposition of the two stationary solutions 
for a = 2.2 and g = 0.5 and different strengths of the loss and gain parameter: 
(a) 7 = 0, (b) 7 = 0.2, (c) 7 = 0.275, (d) 7 = 0.3. 

for that value of g . 

3. Conclusions and outlook 

We have analyzed the simple quantum mechanical model of a Bose- Einstein condensate 
in VT symmetric delta-function double traps by solving the stationary Gross- 
Pitaevskii equation first for negligible nonlinearity and then for finite nonlinearity. 
For vanishing nonlinearity we find that the eigenvalues and eigenstates exhibit the 
same behaviour as in an analogous optical system of two VT symmetric coupled wave 
guides studied by Klaiman et al. [4]. There exist two eigenvalues which merge at a 
critical value of the loss/gain parameter in a branch point, up to the branch point 
the corresponding two eigenstates are VT symmetric, they coalesce at the branch 
point, beyond the branch point the eigenvalues become complex conjugate, and the 
VT symmetry of the wave functions is broken. 

For finite nonlinearity we also find two branches of real eigenvalues which merge 
in an exceptional point. We have the important result that the wave functions of 
the stationary states remain VT symmetric even if the nonlinearity is present. As 
a consequence their moduli are even functions, and therefore the nonlinear state- 
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dependent Hamiltonian selects as solutions exactly such states which make itself VT 
symmetric. We also find a branch of two complex conjugate eigenvalues for which 
the VT symmetry of the wave functions is broken, and with it that of the nonlinear 
Hamiltonian. 

A surprising result is that, with the nonlinearity present, the branches of complex 
conjugate eigenvalues do not bifurcate from the exceptional point, but emerge at a 
smaller value of the loss/gain parameter from the eigenvalue branch of the ground 
state. This implies that there exist ranges of the loss/gain parameter where four 
eigenstates exist. This change in the branching scheme of a VT symmetric nonlinear 
system is in agreement with the results found by Graefe et al. [7J in their investigations 
of a Bose-Hubbard dimcr. 

The puzzle that at the exceptional point no branches of complex eigenvalues 
are born could be solved by analytically continuing the nonlinear term in the Gross- 
Pitaevskii equation, and we could show that indeed two more complex conjugate 
solutions are born at this point. 

The temporal evolution of an initial superposition of the two stationary states 
shows, even for nonvanishing nonlinearity, the typical beating patterns known from 
VT symmetric systems, but we have found that the breakdown of VT symmetry 
occurs at values of the loss/gain parameter smaller than those values where the two 
branches of complex conjugate eigenvalues bifurcate from the ground state or where 
the exceptional point appears. This is due to the time variation of the nonlinear term 
in the Gross-Pitaevskii equation. 

We have considered the case of an attractive nonlinearity but found that the same 
behaviour occurs for repulsive nonlinearity. 

Our investigations can be extended in several directions. First it will be a 
worthwhile enterprise to investigate VT symmetric Bose-Einstein condensates in 
realistic double well potentials [26], in one or more dimensions, and possibly pin down 
physical parameters where the breaking of VT symmetry could be observed in a 
real experiment. Also, in addition to the nonlinearity resulting from the short-range 
contact interaction, condensates with a long-range dipole-dipole interaction [37] could 
be considered. It would also be interesting whether for the nonlinear VT symmetric 
Hamiltonian considered in the present paper similar to the work of Jones [12] the 
construction of a metric operator is possible with respect to which the nonlinear 
Hamiltonian is quasi-Hermitian. Furthermore it would be desirable to find simple 
matrix models which show the unusual branching scheme of the eigenvalues found for 
finite nonlinearity. 
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Appendix A. Gross-Pitaevskii equation 

Expressed in terms of physical quantities the Gross-Pitaevskii equation for a 
condensate in a VT symmetric delta-functions double well reads 
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-G|*(x)| 2 ]*(x) =fx^(x). (A.l) 

The depths of the double wells is determined by Vq , the size of the gain and loss terms 
by ±r. The amplitude of the nonlinear term which arises from the contact interaction 
between the condensate atoms is given by the quantity G and is proportional to the 
scattering length of the condensate atoms, fi is the chemical potential. Measuring 
lengths in units of L = h 2 /(2mVo) and energies in units of E = 2mV" 2 /fi. 2 , brings 
(|A.1|) into the dimensionless form given in ([1]). 

References 

[1] A. Guo, G.J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G.A. 
Siviloglou, and D. N. Christodoulides. Observation of PT-symmetry breaking in complex 
optical potentials. Phys. Rev. Lett, 103:93902, 2009. 
[2] C. E. Riiter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip. 

Observation of parity-time symmetry in optics. Nat. Phys., 6:192, 2010. 
[3] C. M. Bender and S. Boettcher. Real spectra in non-Hermitian Hamiltonians having PT 

symmetry. Phys. Rev. Lett, 80(24) :5243-5246, 1998. 
[4] S. Klaiman, U. Giinther, and N. Moiseyev. Visualization of branch points in PT-symmetric 

waveguides. Phys. Rev. Lett, 101:080402, 2008. 
[5] E.-M. Graefe, H. J. Korsch, and A. E. Niederle. Mean-field dynamics of a non-Hermitian Bose- 

Hubbard dimer. Phys. Rev. Lett, 101:150408, 2008. 
[6] E.-M. Graefe, Uwe Giinther, H. J. Korsch, and A. E. Niederle. A non-Hermitian PT symmetric 
Bose-Hubbard model: eigenvalue rings from unfoldings higher-order exceptional points. J. 
Phys. A, 41:255206, 2008. 
[7] E.-M. Graefe, H. J. Korsch, and A. E. Niederle. Quantum-classical correspondence for a non- 
Hermitian Bose-Hubbard dimer. Phys. Rev. A, 82:013629, 2010. 
[8] H. Ramczani, T. Kottos, R. El-Ganainy, and D. N. Christodoulides. Unidirectional nonlinear 

PT-symmetric optical structures. Phys. Rev. A, 82:043803, 2010. 
[9] K. Rapedius and H. J. Korsch. Resonance solutions of the nonlinear Schroedinger equation in 
an open double-well potential. J. Phys. B, 42:044004, 2009. 
[10] G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher. Symmetry breaking 

in symmetric and asymmetric double-well potentials. Phys. Rev. E, 74:056608, Nov 2006. 
[11] D. Witthaut, K. Rapedius, and H. J. Korsch. The nonlinear Schroedinger equation for the 
delta-comb potential: quasi-classical chaos and bifurcations of periodic stationary solutions. 
J. Nonlin. Math. Phys., 16:207, 2009. 
[12] H. F. Jones. Interface between Hermitian and non-Hermitian Hamiltonians in a model 

calculation. Phys. Rev. D, 78:065032, 2008. 
[13] V. Jakubsky and M. Znojil. An explicitly solvable model of the spontaneous PT-symmetry 

breaking. Czech. J. Phys., 55(9):1113-1116, 2005. 
[14] F. Kh. Abdullaev, A. Abdumalikov, and R. Galimzyanov. Gap solitons in Bose-Einstein 

condensates in linear and nonlinear optical lattices. Phys. Lett. A, 367:149, 2007. 
[15] F. Kh. Abdullaev, A. Gammal, H. L. F. da Luz, and L. Tomio. Dissipative dynamics of matter- 
wave solitons in a nonlinear optical lattice. Phys. Rev. A, 76:43611, 2007. 
[16] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani. Beam dynamics in 

PT symmetric optical lattices. Phys. Rev. Lett, 100:103904, 2008. 
[17] Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides. Optical solitons in 

PT periodic potentials. Phys. Rev. Lett, 100:30402, 2008. 
[18] Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides. Analytical solutions 
to a class of nonlinear Schroedinger equations with PT-like potentials. J. Phys. A, 41:244019, 
2008. 

[19] F. Kh. Abdullaev, V. V. Konotop, M. Salerno, and A. V. Yulin. Dissipative periodic waves, 
solitons, and breathers of the nonlinear Schroedinger equation with complex potentials. Phys. 
Rev. E, 82:56606, 2010. 

[20] Yu. V. Bludov and V. V. Konotop. Nonlinear patterns in Bose-Einstein condensates in 
dissipative optical lattices. Phys. Rev. A, 81:13625, 2010. 

[21] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani. PT-symmetric 
periodic optical potentials. Int. J. Theor. Phys., 50:1019, 2011. 

[22] F. Kh. Abdullaev, Y. V. Kartashov, V. V. Konotop, and D. A. Zezyulin. Solitons in PT- 
symmetric nonlinear lattices. Phys. Rev. A, 83:41805, 2011. 



Nonlinear Schrddinger equation for a VT symmetric delta-functions double well 16 



[23] L. P. Pitaevskii and S. Stringari. Bose-Einstein Condensation. Oxford University Press, 2003. 
[24] L. D. Carr, C. W. Clark, and W. P. Reinhardt. Stationary solutions of the one-dimensional 

nonlinear Schroedingcr equation. I. Case of repulsive nonlinearity. Phys. Rev. A, 62:63610, 

2000. 

[25] L. D. Carr, J. N. Kutz, and W. P. Reinhardt. Stability of stationary states in the cubic 
nonlinear Schroedingcr equation: Applications to the Bosc-Einstcin condensate. Phys. Rev. 
E, 63:66604, 2001. 

[26] D. Ananikian and T. Bcrgcman. Gross-Pitaevskii equation for Bose particles in a double-well 
potential: Two-mode models and beyond. Phys. Rev. A, 73:13604, 2006. 

[27] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau. The physics of dipolar bosonic 
quantum gases. Rep. Progr. Phys., 72:126401, 2009. 



